#!/bin/tcsh -ef
####
#
#############################################################################
#                                                                           #
# Title: Get Lattice & Tilt                                                 #
#                                                                           #
# (C) 2dx.org, GNU Plublic License.                                         #
#                                                                           #
# Created..........: 02/20/2006                                             #
# Last Modification: 09/05/2007                                             #
# Author...........: 2dx.org                                                #
#                                                                           #
#############################################################################
#
# SORTORDER: 30
#
# MANUAL: This script determines the <B>reciprocal crystal lattice</B> in a Fourier transformation of the image. Two different algorithms are offered. The first algorithm <I>FindLattice</I> makes use of a-prior knowledge of the real-space dimensions of the crystal lattice, and of the tilt geometry. The second algorithm <I>GetLattice</I> performs a peak list difference calculation and determines the most likely lattice from the power-spectrum.
# 
# MANUAL: <B><U>FindLattice</U></B>
#
# MANUAL: This algorithm determines the lattice, using the a-priori knowledge of the expected lattice and the approximate tilt geometry (from the defocus determination). Once the lattice is determined, its distortion is used to determine the tilt geometry again, using the MRC program EMTILT (D. Agard). With tilt angles higher than 25 degrees, the lattice distortion will usually allow a better determination of the tilt geometry than from a defocus gradient. 
#
# MANUAL: <B><U>GetLattice</U></B>
#
# MANUAL: This algorithm <I>Difference Lattice</I> does not require a priori knowledge of the real-space lattice. 
#
# MANUAL: <i>2dx_diffLattice</i> takes advantage of the shifted origin peak list. The results of this shifting are identical to those produced by weighted, difference-vector methods --with the added benefit of using an adaptive criteria for peak-merging as opposed to standard Gaussian determination.
#
# MANUAL: <i>2dx_diffLattice</i> begins by finding all peaks with radii less than the maximum radius of the 8 strongest peaks from the list. It follows from the basis for the shifted origin that the strongest peak will necessarily fall on the stronger of the two lattices. Every possible vector pair is then used to generate a lattice which is then compared against the available peak list. Lattice error is computed by transforming peaks into Miller indices, then measuring their divergence from integer values, producing a result which is then multiplied by the strength of the point in question. For lattice pairs with an angle of 10 degrees or less between them, a further weighting factor of Exp(10-Theta) is multiplied, with theta the measured angle between the vectors in degrees. (This somewhat arbitrary weighting is used to prevent extremely high density lattices from overpowering sparse, yet more accurate, ones.)
#
# MANUAL: The vector pair which generates the lowest error lattice is then used as a Miller basis. Using this basis, a Miller index is then assigned to any peak which lies within epsilon=0.15*sqrt(2) of a point on the lattice. A peak assigned to an index is replaced if another peak is found to lie closer to the lattice point. Peaks which are not assigned lattice points are discarded from further consideration.
#
# MANUAL: Using the generated Miller index - peak position pairs, a least squares fit is then performed to determine a candidate lattice which is then used as the basis for iterative refinement, according to the above procedure, to produce the <B>final lattice</B>.
#
# MANUAL: <HR>
#
# MANUAL: The here determined tilt geometry will be entered in the status pane in the second column, under "Latt.". 
#
# MANUAL: This lattice-based tilt geometry is only reliable for larger tilt angles, which should be greater than 25 degrees. For smaller tilt angles, the defocus gradient is the safer method. In that case, the script "Get Sample Tilt" will use the now determined lattice orientation to translate the defocus-based tilt geometry into the sample tilt geometry.
#
#
# PUBLICATION: 2dx - User-friendly image processing for 2D crystals: <A HREF="http://dx.doi.org/10.1016/j.jsb.2006.07.020">J. Struct. Biol. 157 (2007) 64-72</A>
# PUBLICATION: Image processing of 2D crystal images: <A HREF="http://link.springer.com/protocol/10.1007%2F978-1-62703-176-9_10">Methods in Molecular Biology (2012)</A>
#
# DISPLAY: imagesidelength
# DISPLAY: TLTAXIS
# DISPLAY: TLTANG
# DISPLAY: TLTAXA
# DISPLAY: TAXA
# DISPLAY: TANGL
# DISPLAY: lattice
# DISPLAY: secondlattice
# DISPLAY: secondlattice_exists
# DISPLAY: secondlattice_clone
# DISPLAY: secondlattice_is_this
# DISPLAY: backuplattice
# DISPLAY: defocus
# DISPLAY: beamtilt
# DISPLAY: defocusfield
# DISPLAY: RESMAX
# DISPLAY: realang
# DISPLAY: realcell
# DISPLAY: SYM
# DISPLAY: tempkeep
# DISPLAY: delta
# DISPLAY: streakfactor
# DISPLAY: FirstPeakNum
# DISPLAY: peakNum
# DISPLAY: testHand
# DISPLAY: regenPL
# DISPLAY: InnerExclusionRadius
# DISPLAY: comment
# DISPLAY: det_lattice
# DISPLAY: det_tilt_lat
# DISPLAY: do_lattice_algorithm
# DISPLAY: laterror_abserror
# DISPLAY: laterror_rmsderror
# DISPLAY: laterror_peaksused
# DISPLAY: laterror_nodepeakdensity
#
#$end_local_vars
#
#
set bin_2dx = ""
set proc_2dx = ""
#
set SYN_Unbending = ""
set PHASEORI_done = ""
#
set imagename = ""
set nonmaskimagename = ""
set imagenumber = ""
set realcell = ""
set lattice = ""
set TLTAXIS = ""
set TLTANG = ""
set TLTAXA = ""
set TAXA = ""
set TANGL = ""
set imagesidelength = ""
set secondlattice = ""
set secondlattice_exists = ""
set secondlattice_is_this = ""
set secondlattice_clone = ""
set magnification = ""
set stepdigitizer = ""
set tempkeep = ""
set revhk = ""
set realang = ""
set RESMIN = ""
set RESMAX = ""
set ALAT = "" 
set det_tilt_lat = ""
set delta = ""
set streakfactor = ""
set FirstPeakNum = ""
set peakNum = ""
set testHand = ""
set regenPL = ""
set InnerExclusionRadius = ""
set det_lattice = ""
set do_lattice_algorithm = ""
set laterror_abserror = ""
set laterror_rmsderror = ""
set laterror_peaksused = ""
set laterror_nodepeakdensity = ""
set MASKING_done = ""
set LATTICE_done = ""
set defocus = ""
#
#$end_vars
#
set scriptname = 2dx_getLattice
#
\rm -f LOGS/${scriptname}.results
#
set IS_2DX = yes
source ${proc_2dx}/initialize
#
echo "<<@evaluate>>"
#
set docfile = 'peaks_xy_final.dat'
set docfileraw = 'peaks_xy.dat'
set iradius = '8'
set imaxrad = '20'
#
set date = `date`
echo date = ${date}
#
if ( ${det_lattice}x == 'x' ) then
  set det_lattice = 'y'
  echo "set det_lattice = ${det_lattice}" >> LOGS/${scriptname}.results
  echo ":: det_lattice corrected to ${det_lattice}"
endif
#
if ( x${TLTAXIS} == "x-" || x${TLTAXIS} == "x--" ) then
  set TLTAXIS = "0.0"
endif
if ( x${TLTANG} == "x-" || x${TLTANG} == "x--" ) then
  set TLTANG = "0.0"
endif
#
echo "::Defocus = "${defocus}
echo "# IMAGE-IMPORTANT: "FFTIR/${imagename}_fft.mrc "<FFT of Image>" >> LOGS/${scriptname}.results
echo "# IMAGE-IMPORTANT: "FFTIR/${imagename}_red_fft.mrc "<FFT of Downsampled Image>" >> LOGS/${scriptname}.results
#
if ( ${det_lattice} == 'n' && ${det_tilt_lat} == 'n' ) then
  ${proc_2dx}/linblock "Skipping. Not determining lattice."
  exit
endif
#
if ( (${SYN_Unbending} != "0") && (${PHASEORI_done} == "y") ) then
  ${proc_2dx}/linblock "Skipping. Synthetical Reference is used."
  ${proc_2dx}/linblock "Lattice should be determined already."
  echo "# IMAGE: SCRATCH/${imagename}_taper.mrc <Edge-Tapered Image>" >> LOGS/${scriptname}.results
  echo "# IMAGE: SCRATCH/pks-non-masked_image.mrc <Powerspectrum>" >> LOGS/${scriptname}.results
  echo "# IMAGE: SCRATCH/pks-amp_LHpass.mrc <Band-Passed Powerspectrum>" >> LOGS/${scriptname}.results
  echo "# IMAGE: SCRATCH/pks-masked_image.mrc <Masked and Band-Passed Powerspectrum>" >> LOGS/${scriptname}.results
  echo "# IMAGE-IMPORTANT: SCRATCH/pks-average.mrc <Origin-Shifted Average Powerspectrum>"  >> LOGS/${scriptname}.results
  exit
endif
#
\rm -f TMP456789.tmp
#
echo "<<@progress: 1>>"
#
if ( ${det_lattice} == 'y' ) then
  #
  set ialg = `echo ${do_lattice_algorithm} | cut -c1-1` 
  #
  if ( ${ialg} == "0" ) then
    #
    #############################################################################
      ${proc_2dx}/linblock "Running FindLattice"
    #############################################################################
    #
    if ( "${realcell}" == "100.0,100.0" || "${realcell}" == "100,100" ) then
      echo ":: "
      echo ":: "
      echo "::                        WARNING: You are using FindLattice,"
      echo "::                     but did you define a real lattice already?" 
      echo ":: "
      echo "::        (The current value of ${realcell} looks like you didn't do that yet.)"
      echo "::       You might need to first determine the real unit cell length and angle,"
      echo "::                            or use GetLattice instead."
      echo ":: "
      echo ":: "
      echo "#WARNING: Warning: For FindLat you need a valid real lattice first." >> LOGS/${scriptname}.results
    endif
    #
    if ( "${realcell}" == "0,0" || "${realcell}" == "0.0,0.0" ) then
      echo ":: "
      echo ":: "
      echo "::                        ERROR: You are using FindLattice,"
      echo "::                     but there is no real lattice defined." 
      echo ":: "
      echo "::            (The current value for the real lattice is ${realcell}.)"
      echo "::        You need to first determine the real unit cell length and angle."
      echo "::      You can do that manually, or use GetLattice instead. After that, run"
      echo "::     the special script Evaluate Lattice to determine the real cell values."
      echo ":: "
      echo ":: "
      echo "#WARNING: ERROR: For FindLat you need a valid real lattice first." >> LOGS/${scriptname}.results
      ${proc_2dx}/protest "Aborting."
    endif
    #
  else
    #
    #############################################################################
    ${proc_2dx}/linblock "Running GetLattice"
    #############################################################################
    #
  endif
  echo "<<@progress: 5>>"
  #
  if ( ${regenPL} == "y" || ( ! -e ${docfile} ) || ( ! -e ${docfileraw} ) ) then
    #
    #############################################################################
    ${proc_2dx}/linblock "2dx_taperedge: To taper the edges to prevent stripes in the FFT"
    #############################################################################
    #
    set largeimage = `echo ${imagesidelength} | awk '{ if ( $1 > 6000 ) { s = 1 } else { s = 0 } } END { print s }'`
    #
    if ( ${largeimage} == '1' ) then
      if ( ! -e FFTIR/${imagename}_red.mrc ) then
        ${proc_2dx}/protest "::ERROR: Downsampled image does not exist. First run Calculate-FFT."
      endif
      setenv IN FFTIR/${imagename}_red.mrc
      ${proc_2dx}/lin "Doing peaksearch on downsampled image"
    else
      setenv IN ${imagename}.mrc
      ${proc_2dx}/lin "Doing peaksearch on full-sized image"
    endif
    # 
    \rm -f     SCRATCH/${imagename}_taper.mrc
    setenv OUT SCRATCH/${imagename}_taper.mrc
    ${bin_2dx}/2dx_taperedgek.exe << eot
10,10,100,10       ! IAVER,ISMOOTH,ITAPER 
eot
    #
    echo "<<@progress: 10>>"
    #
    #############################################################################
    ${proc_2dx}/linblock "2dx_peaksearch: To generate a completed peak list ..."
    #############################################################################
    #
    echo "<<@progress: 20>>"
    #
    \rm -f ${docfile}
    \rm -f ${docfileraw}
    \rm -f 2dx_peaksearch-amp_before_LHpass.mrc
    \rm -f 2dx_peaksearch-amp_LHpass.mrc
    #
    echo "FirstPeakNum = ${FirstPeakNum}"
    echo "peakNum = ${peakNum}"
    echo "InnerExclusionRadius = ${InnerExclusionRadius}"
    ${bin_2dx}/2dx_peaksearch.exe SCRATCH/${imagename}_taper.mrc ${FirstPeakNum} ${peakNum} ${InnerExclusionRadius} ${imagesidelength} ${streakfactor}
    #endif
    #
    if ( ${tempkeep} == 'y' ) then
      #
      \mv -f 2dx_peaksearch-non-masked_image.mrc SCRATCH/pks-non-masked_image.mrc
      if ( -e 2dx_peaksearch-amp_before_LHpass.mrc ) then
        \mv -f 2dx_peaksearch-amp_before_LHpass.mrc SCRATCH/pks-amp_before_LHpass.mrc
      endif
      if ( -e 2dx_peaksearch-amp_LHpass.mrc ) then
        \mv -f 2dx_peaksearch-amp_LHpass.mrc SCRATCH/pks-amp_LHpass.mrc
      endif
      \mv -f 2dx_peaksearch-masked_image.mrc SCRATCH/pks-masked_image.mrc
      \rm -f 2dx_peaksearch-phase.mrc
      echo "# IMAGE: SCRATCH/${imagename}_taper.mrc <Edge-Tapered Image>" >> LOGS/${scriptname}.results
      echo "# IMAGE: SCRATCH/pks-non-masked_image.mrc <Powerspectrum>" >> LOGS/${scriptname}.results
      echo "# IMAGE: SCRATCH/pks-amp_before_LHpass.mrc <Pre Band-Pass Powerspectrum>" >> LOGS/${scriptname}.results
      echo "# IMAGE: SCRATCH/pks-amp_LHpass.mrc <Post Band-Pass Powerspectrum>" >> LOGS/${scriptname}.results
      echo "# IMAGE: SCRATCH/pks-masked_image.mrc <Masked and Band-Passed Powerspectrum>" >> LOGS/${scriptname}.results
      #
    else
      \rm -f SCRATCH/${imagename}_taper.mrc
      \rm -f 2dx_peaksearch-non-masked_image.mrc
      \rm -f 2dx_peaksearch-amp_LHpass.mrc
      \rm -f 2dx_peaksearch-masked_image.mrc
      \rm -f 2dx_peaksearch-phase.mrc
    endif
    #
    \mv -f 2dx_peaksearch-average.mrc SCRATCH/pks-average.mrc
    echo "# IMAGE-IMPORTANT: SCRATCH/pks-average.mrc <Origin-Shifted Average Powerspectrum>"  >> LOGS/${scriptname}.results
    #
  else
    #############################################################################
    ${proc_2dx}/linblock "2dx_peaksearch: Skipping"
    #############################################################################
    echo "# IMAGE: SCRATCH/${imagename}_taper.mrc <Edge-Tapered Image>" >> LOGS/${scriptname}.results
    echo "# IMAGE: SCRATCH/pks-non-masked_image.mrc <Powerspectrum>" >> LOGS/${scriptname}.results
    echo "# IMAGE: SCRATCH/pks-amp_LHpass.mrc <Band-Passed Powerspectrum>" >> LOGS/${scriptname}.results
    echo "# IMAGE: SCRATCH/pks-masked_image.mrc <Masked and Band-Passed Powerspectrum>" >> LOGS/${scriptname}.results
    echo "# IMAGE-IMPORTANT: SCRATCH/pks-average.mrc <Origin-Shifted Average Powerspectrum>" >> LOGS/${scriptname}.results
    #
  endif
  #
  echo "<<@evaluate>>"
  #
  if ( ${ialg} == "0" ) then
    #
    echo "<<@progress: 35>>"
    #
    #############################################################################
    ${proc_2dx}/linblock "2dx_findlat: To determine a lattice ..."
    #############################################################################
    #
    \rm -f 2dx_findlat.out
    # 
    set tltAx = ${TLTAXIS}
    set tltAn = ${TLTANG}
    #
    if ( ${testHand} == 'n' ) then
      set itesthand = 0
    else
      set itesthand = 1
    endif
    #  
    echo "starting findlat with:"
    echo "${realcell}"
    echo "${realang}"
    echo "${tltAx}"
    echo "${tltAn}"
    echo "${imagesidelength}"
    echo "${magnification}"
    echo "${stepdigitizer}"
    echo "${docfile}"
    echo "0.1,0.008,8,${delta},16,300,${itesthand}"
    echo " "
    #
  ${bin_2dx}/2dx_findlat.exe << eot
${realcell}
${realang}
${tltAx}
${tltAn}
${imagesidelength}
${magnification}
${stepdigitizer}
${docfile}
0.1,0.008,8,${delta},16,300,${itesthand} !AngleStepSize,MagVariation,Nmin,Delta,MaxHK,MaxSpot
eot
    #
    echo "findlat finished."
    #
    if ( ${tempkeep} == 'n' ) then
      #\rm -f peaks_xy_final.dat
    endif
    #
    if ( ! -e 2dx_findlat.out ) then
      ${proc_2dx}/protest "2dx_findlat: ERROR occured."
    endif
    #
    set tmplat = `cat 2dx_findlat.out | head -n 1`
    set firstscore  = `cat 2dx_findlat.out | head -n 2 | tail -n 1`
    set newlat = `echo ${tmplat} | sed 's/ /,/g'`
    set oldlat = `echo ${lattice}`
    #
    echo "set lattice = "\"${newlat}\" >> LOGS/${scriptname}.results
    set lattice = ${newlat}
    echo "::  Old lattice is ${oldlat}"
    echo "::  New lattice is ${lattice} (Score = ${firstscore})"
    #
    echo " " >> History.dat
    echo ":Date: ${date}" >> History.dat
    echo "::Lattice is ${lattice}" >> History.dat
    #
    if ( ${secondlattice_is_this} == "n" ) then
      set tmplat = `cat 2dx_findlat.out | head -n 3 | tail -n 1`
      set secondscore  = `cat 2dx_findlat.out | head -n 4 | tail -n 1`
      set newsecondlattice = `echo ${tmplat} | sed 's/ /,/g'`
      set oldsecondlattice = `echo ${secondlattice}`
      #
      # \rm -f 2dx_findlat.out
      #
      echo "::  Old 2nd lattice is ${oldsecondlattice}"
      set goodsecond = `echo ${secondscore} | awk '{ if ( $1 > 0.6 ) { s = 1 } else { s = 0 }} END { print s }'`
      if ( ${goodsecond} == "1" ) then
        set secondlattice = ${newsecondlattice}
        echo "::  New 2nd lattice is ${secondlattice} (Score = ${secondscore})"
        echo "set secondlattice = "\"${newsecondlattice}\" >> LOGS/${scriptname}.results
        echo "set secondlattice_exists = y" >> LOGS/${scriptname}.results
        #
        if ( ${secondlattice_clone} == "y" ) then
          #############################################################################
          echo ":: "
          ${proc_2dx}/linblock "Cloning this image directory to initiate processing of second lattice"
          echo ":: "
          #############################################################################
          #
          set secondlattice_clone = "n"
          echo "set secondlattice_clone = ${secondlattice_clone}" >> LOGS/${scriptname}.results
          #
          set olddir = `pwd`
          set currentdir = `echo ${olddir} | rev | cut -d\/ -f1 | rev`
          set newdir = `echo ${currentdir}_lat2`
          cd ..
          if ( ! -d ${newdir} ) then
            \mkdir ${newdir}
            cd ${newdir}
            set new_created_dir = $PWD
            if ( -e  ../${currentdir}/${nonmaskimagename}.tif ) then
              \cp -f ../${currentdir}/${nonmaskimagename}.tif .
            endif
            if ( -e  ../${currentdir}/${nonmaskimagename}.mrc ) then
              \cp -f ../${currentdir}/${nonmaskimagename}.mrc .
            endif
            if ( -e  ../${currentdir}/${nonmaskimagename}_stack.mrc ) then
              \cp -f ../${currentdir}/${nonmaskimagename}_stack.mrc .
            endif
            if ( -e  ../${currentdir}/${nonmaskimagename}_raw.mrc ) then
              \cp -f ../${currentdir}/${nonmaskimagename}_raw.mrc .
            endif
            \cp -f ../${currentdir}/2dx_image.cfg .
            set newimagenumber = `echo ${imagenumber} | awk '{ s = $1 + 1 } END { print s }'`
            echo "set imagenumber = ${newimagenumber}" >> 2dx_image.cfg
            echo "set det_lattice = n" >> 2dx_image.cfg
            echo "set lattice = "\"${newsecondlattice}\" >> 2dx_image.cfg
            echo "set secondlattice = "\"${lattice}\" >> 2dx_image.cfg
            echo "set secondlattice_is_this = "\"y\" >> 2dx_image.cfg
            echo "set secondlattice_exists = "\"y\" >> 2dx_image.cfg
            ${proc_2dx}/linblock "Image cloned, new directory name is ${newdir}"
            #
            cd ${olddir}
          endif
        endif
        #
      else
        set secondlattice = "0.0,0.0,0.0,0.0"
        echo "::  New 2nd lattice would have been ${newsecondlattice} (Score = ${secondscore}),"
        echo "::  Due to the bad score it will not be used."
      endif
    else
      echo "set secondlattice_exists = y" >> LOGS/${scriptname}.results
    endif
    #
  else
    #
    #############################################################################
    ${proc_2dx}/linblock "2dx_getlat: To determine a lattice ..."
    #############################################################################
    #
    set tltAx = ${TLTAXIS}
    set tltAn = ${TLTANG}
    #
    set lattice=`${bin_2dx}/2dx_getlat.exe ${docfile} | tail -n 1`
    #
    echo set lattice = \"${lattice}\" >> LOGS/${scriptname}.results
    #
    echo " " >> History.dat
    echo ":Date: ${date}" >> History.dat
    echo "::Lattice is ${lattice}" >> History.dat
    #
    if ( ${tempkeep} == 'n' ) then
      #\rm -f peaks_xy_final.dat
    endif
    #
    echo "::  New lattice is ${lattice}"
  endif
  #
  echo "<<@progress: 50>>"
  echo "<<@evaluate>>"
  #
  #############################################################################
  ${proc_2dx}/linblock "2dx_getspot - To generate a first spotlist."
  #############################################################################
  #
  ${proc_2dx}/lin "2dx_getspot.exe - peaks_xy.dat + lattice => 2dx_getspot.out"
  #
  ${bin_2dx}/2dx_getspot.exe << eot
${lattice}
${docfileraw}
${iradius}
${imaxrad}
eot
  #
  ${proc_2dx}/lin "cp 2dx_getspot.out ${nonmaskimagename}.spt"
  \cp -f 2dx_getspot.out ${nonmaskimagename}.spt
  \rm -f 2dx_getspot.out
  #
    #############################################################################
    #############################################################################
    #############################################################################
    #############################################################################
else
  if ( "${lattice}" == "0.0,100.0,100.0,0.0" || "${lattice}" == "0.0,0.0,0.0,0.0" ) then
    ${proc_2dx}/protest "ERROR: you need to determine the reciprocal lattice first."
  endif
  if ( "${realcell}" == "100.0,100.0" || "${realcell}" == "100,100" ) then
    echo "::        (The current value of ${realcell} looks like you didn't do that yet.)"
    echo "::       You might need to first determine the real unit cell length and angle,"
    ${proc_2dx}/protest "ERROR: Aborting."
  endif
  if ( "${realcell}" == "0,0" || "${realcell}" == "0.0,0.0" ) then
    echo "::            (The current value for the real lattice is ${realcell}.)"
    echo "::        You need to first determine the real unit cell length and angle."
    echo "::      You can do that manually, or use GetLattice instead. After that, run"
    echo "::     the special script Evaluate Lattice to determine the real cell values."
    ${proc_2dx}/protest "Aborting."
  endif
endif
#############################################################################
#############################################################################
#############################################################################
#############################################################################
#
#############################################################################
${proc_2dx}/linblock "2dx_lencalc: To calculate some values for the tilt geometry determination."
#############################################################################
#
set lencalc_docfile = "2dx_lencalc-doc.tmp"
#
\rm -f ${lencalc_docfile}
#
${bin_2dx}/2dx_lencalc.exe << eot
${lattice}
${realcell}
${realang}
${imagesidelength}
${magnification}
${stepdigitizer}
${lencalc_docfile}
eot
#
cat ${lencalc_docfile}
#
set hand     = `head -n 2 ${lencalc_docfile} | tail -n 1`
set untilted = `head -n 4 ${lencalc_docfile} | tail -n 1`
set alpha    = `head -n 6 ${lencalc_docfile} | tail -n 1`
set tilted   = `head -n 8 ${lencalc_docfile} | tail -n 1`
#
\rm -f ${lencalc_docfile}
#
echo "<<@progress: 60>>"
#
if ( ${det_tilt_lat} == 'y' ) then
  #
  #############################################################################
  #                                                                           #
  #  EMTILT                                                                   #
  #                                                                           #
  #                 ->  1   , 1   , 90.0       (or also : lenA,lenB,90.0)     #
  #                 ->  lenH, lenK, ang(H,K)                                  #
  #                                                                           #
  #                 <-  ang(Tiltaxis to A*) = Omega                           #
  #                 <-  Tiltangle (without sign)                              #
  #                                                                           #
  #  A* was in the EM in the sample-plane.                                    #
  #  Omega is NOT the angle from the Tiltaxis to H, but this one without      #
  #     perspective. Omega is in the negative plane.                          #
  #     Omega is slightly smaller than the latter (TLTAXA).                   #
  #                                                                           #
  #############################################################################
  #                                                                           #
  #  Angle from TLT-axis to H is TLTAXA, TLT-axis to K is TLTAXB. (???)       #
  #  Angle from Y-axis to H is TLTAXA, Y-axis to K is TLTAXB.                 #
  #  Henderson: "TLTANG, TLTAXA, TLTAXB are on negative, nobody wants them."  #
  #  This is not fully correct: TTMASK or TTBOX need them.                    #
  #                                                                           #
  #  TLTAXIS: Angle from X-axis to Tiltaxis                                   #
  #  TLTANG is positive for strong underfocus at top of image.                #
  #  TLTANG is positive for left overfocus, right underfocus when TLTAXIS= Y  #
  #  TLTANG is positive for strong underfocus on left side, when TLTAXIS= Y   #
  #       (Left screw positive counting)                                      #
  #          TLTAXIS, TLTANG are for TTREFINE.                                #
  #                                                                           #
  #  TTREFINE has perfect output for ORIGTILTC :                              #
  #          TAXA  , TAXB  , TANGL  are on sample, for ORIGTILTC.             #
  #                                                                           #
  #############################################################################
  #
  # Henderson in an email 6/2000:
  #   y=0 means the x-axis which is horizontal.  It is indeed the bottom line of
  # the image.  The first point on that line is (0,0), the second point (1,0)
  # and so on [e.g. (x,0)].  Therefore all points on the line have y=0.  If the
  # tiltaxis is horizontal, parallel to x and the defocus gets stronger as you
  # go to higher values of y (i.e. further up the image), then TLTANGL is positive.
  # This is a robust definition and TLTAXIS can take any number between -89.999
  # and +89.999.  If it happens to be exactly 90.000 (note that this has never
  # occurred in practice), then you have to know how the program reacts and that is  
  # where the extra x=0 description arises.  In practice, if your film orientation
  # is such that the tiltaxis is mostly vertical but varies a bit with some of
  # them on different sides of the vertical, then those with TLTAXIS of 85 degrees,
  # for example, would have TLTANG positive, whereas those with TLTAXIS of 95 (i.e.
  # -85) degrees would have TLTANGL negative.  So, the TLTANGL sign would change
  # at 90.   This does not happen with the tiltaxis roughly horizontal, when there
  # is no change in the definition as the TLTAXIS passes from -10 to +10.  I hope
  #?| this is now clearer.
  #
  #############################################################################
  #
  #############################################################################
  ${proc_2dx}/linblock "emtilt: To calculate the tilt geometry from lattice distortions."
  #############################################################################
  #
  \rm -f TMP.2dx_emtilt.1.out
  #
  if ( ${TLTAXIS} == "-" || ${TLTAXIS} == "--" ) then
    set oldTLTAXIS = "0.0"
  else
    set oldTLTAXIS = ${TLTAXIS}
  endif
  if ( ${TLTANG} == "-" || ${TLTANG} == "--" ) then
    set oldTLTANG = "0.0"
  else
    set oldTLTANG  = ${TLTANG}
  endif
  #
  set emtilt_docfile = "2dx_emtilt-doc.tmp"
  \rm -f ${emtilt_docfile} 
  #
  echo untilted = ${untilted}
  echo tilted = ${tilted}
  echo oldTLTAXIS = ${oldTLTAXIS}
  echo oldTLTANG = ${oldTLTANG}
  echo hand = ${hand}
  echo alpha = ${alpha}
  #
  ${bin_2dx}/2dx_emtilt.exe << eot
${untilted}
${tilted}
${oldTLTAXIS}
${oldTLTANG}
${hand}
${alpha}
${emtilt_docfile}
eot
  #
  echo "2dx_emtilt.exe finished."  
  #
  # read those values for the tilt geometry back into this script shell:
  set LATTICE_TLTAXIS = `head -n 2  ${emtilt_docfile} | tail -n 1`
  set LATTICE_TLTANG  = `head -n 4  ${emtilt_docfile} | tail -n 1`
  set LATTICE_TLTAXA  = `head -n 6  ${emtilt_docfile} | tail -n 1`
  set LATTICE_TAXA    = `head -n 8  ${emtilt_docfile} | tail -n 1`
  set LATTICE_TANGL   = `head -n 10 ${emtilt_docfile} | tail -n 1`
  \rm -f ${emtilt_docfile}  
  #
  echo "<<@progress: 75>>"
  #
  set highangle = `echo ${LATTICE_TLTANG} | awk '{if($1>25 || $1<-25) {s=1} else {s=0}} END {print s}'`
  #
else
  set highangle = 0
endif
#
if ( ${det_tilt_lat} == 'y' ) then
  #
  if ( ${highangle} == '1' ) then
    #############################################################################
    ${proc_2dx}/linblock "Saving tilt geometry from lattice, because of high tilt."
    #############################################################################
    set TLTAXIS = ${LATTICE_TLTAXIS}
    set TLTANG  = ${LATTICE_TLTANG}
    set TLTAXA  = ${LATTICE_TLTAXA}
    set TAXA    = ${LATTICE_TAXA}
    set TANGL   = ${LATTICE_TANGL}
    #
    # Make sure these values are going back into 2dx:
    echo "set TLTAXIS = ${TLTAXIS}" >> LOGS/${scriptname}.results
    echo "set TLTANG  = ${TLTANG}"  >> LOGS/${scriptname}.results
    echo "set TLTAXA  = ${TLTAXA}"  >> LOGS/${scriptname}.results
    echo "set TAXA    = ${TAXA}"    >> LOGS/${scriptname}.results
    echo "set TANGL   = ${TANGL}"   >> LOGS/${scriptname}.results
    echo "set DEFOCUS_ACTIVE = 2"   >> LOGS/${scriptname}.results
    #
  else
    #############################################################################
    echo ":: "
    ${proc_2dx}/linblock "NOT saving tilt geometry from lattice, because of low tilt."
    echo ":: "
    #############################################################################
    echo "TLTAXIS = ${LATTICE_TLTAXIS}"
    echo "TLTANG  = ${LATTICE_TLTANG}" 
    echo "TLTAXA  = ${LATTICE_TLTAXA}" 
    echo "TAXA    = ${LATTICE_TAXA}"   
    echo "TANGL   = ${LATTICE_TANGL}"  
    echo "(These values were calculated, but will not be used.)"
  endif
  #
  echo "set LATTICE_TLTAXIS = "\"${LATTICE_TLTAXIS}\" >> LOGS/${scriptname}.results
  echo "set LATTICE_TLTANG  = "\"${LATTICE_TLTANG}\"  >> LOGS/${scriptname}.results
  echo "set LATTICE_TLTAXA  = "\"${LATTICE_TLTAXA}\"  >> LOGS/${scriptname}.results
  echo "set LATTICE_TAXA    = "\"${LATTICE_TAXA}\"    >> LOGS/${scriptname}.results
  echo "set LATTICE_TANGL   = "\"${LATTICE_TANGL}\"   >> LOGS/${scriptname}.results
  #
  echo " " >> History.dat
  echo "::From Lattice:" >> History.dat
  echo "::TLTAXIS = ${LATTICE_TLTAXIS}" >> History.dat
  echo "::TLTANG  = ${LATTICE_TLTANG}" >> History.dat
  echo "::TLTAXA  = ${LATTICE_TLTAXA}" >> History.dat
  echo "::TAXA    = ${LATTICE_TAXA}" >> History.dat
  echo "::TANGL   = ${LATTICE_TANGL}" >> History.dat
  #
endif
#
#############################################################################
${proc_2dx}/linblock "2dx_calcmag - to calculate the theoretical magnification"
#############################################################################
echo "<<@progress: 85>>"
echo "<<@evaluate>>"
#
set outputfile = 2dx_verifyParams.tmp
\rm -f ${outputfile}
#
echo ": running ${bin_2dx}/2dx_calcmag.exe  with:"
echo ":"${realcell}
echo ":"${realang}
echo ":"${TLTAXIS}
echo ":"${TLTANG}
echo ":"${lattice}
echo ":"${imagesidelength}
echo ":"${magnification}
echo ":"${stepdigitizer}
echo ":"${outputfile}
echo ":"
#
${bin_2dx}/2dx_calcmag.exe << eot
${realcell}
${realang}
${TLTAXIS}
${TLTANG}
${lattice}
${imagesidelength}
${magnification}
${stepdigitizer}
${outputfile}
eot
#
if ( ! -e ${outputfile} ) then
  ${proc_2dx}/protest "ERROR in 2dx_calcmag.exe"
endif
#
set theormag = `cat ${outputfile} | head -n 1`
set RANGrec  = `cat ${outputfile} | head -n 2 | tail -n 1`
set RANGreal = `cat ${outputfile} | head -n 3 | tail -n 1`
#
\rm -f ${outputfile}
#
${proc_2dx}/linblock "Theoretical magnification is ${theormag}, given magnification is ${magnification}"
${proc_2dx}/linblock "Angle in reciprocal lattice is ${RANGrec}."
${proc_2dx}/linblock "Angle in real-space lattice is ${RANGreal}."
#
echo "set CALCULATEDMAG = ${theormag}" >> LOGS/${scriptname}.results
#
#############################################################################
${proc_2dx}/linblock "2dx_laterror to find fitness of current lattice."
#############################################################################
#
${bin_2dx}/2dx_laterror.exe "${lattice}" ${docfile} > SCRATCH/2dx_laterror.tmp
cat SCRATCH/2dx_laterror.tmp
set laterror_abserror=`cat SCRATCH/2dx_laterror.tmp | tail -n 7 | head -n 1`
set laterror_rmsderror=`cat SCRATCH/2dx_laterror.tmp | tail -n 5 | head -n 1`
set laterror_peaksused=`cat SCRATCH/2dx_laterror.tmp | tail -n 3 | head -n 1 | cut -f1 -d " "`
set laterror_nodepeakdensity=`cat SCRATCH/2dx_laterror.tmp | tail -n 1`
#
echo ": "
echo "::Absolute Error: ${laterror_abserror}"
echo "::RMSDn ERROR: ${laterror_rmsderror}"
echo "::Peaks Used: ${laterror_peaksused}"
echo "::Nodes/Peak Density: ${laterror_nodepeakdensity}"
#
echo "set laterror_abserror = ${laterror_abserror}" >> LOGS/${scriptname}.results
echo "set laterror_rmsderror = ${laterror_rmsderror}" >> LOGS/${scriptname}.results
echo "set laterror_peaksused = ${laterror_peaksused}" >> LOGS/${scriptname}.results
echo "set laterror_nodepeakdensity = ${laterror_nodepeakdensity}" >> LOGS/${scriptname}.results
#
echo "set LATTICE_done = y" >> LOGS/${scriptname}.results
echo "set SPOTS_done = n" >> LOGS/${scriptname}.results
echo "set UNBENDING_done = n" >> LOGS/${scriptname}.results
echo "set ML_done = n" >> LOGS/${scriptname}.results
echo "set CTF_done = n" >> LOGS/${scriptname}.results
echo "set MERGING_done = n" >> LOGS/${scriptname}.results
#
echo "<<@progress: 100>>"
#
##########################################################################
${proc_2dx}/linblock "${scriptname} - normal end."
##########################################################################
#
